#!/usr/bin/env python
# -*- coding: utf-8 -*-

""" By Martin Senande-Rivera
    For Towards and atmosphere more favourable to firestorm development in Europe """

import os
import glob, sys
import numpy as np
import xarray as xr
import pandas as pd
from scipy import stats

GCM=sys.argv[1]  # Enter GCM
RCM=sys.argv[2]  # Enter RCM
rcp=sys.argv[3]  # Enter RCP scenario

path_outs='./'+RCM+'/'+GCM+'/'

# Number of days with FWI above the threshold
FWI = xr.open_dataset(path_outs+'FWI_ndays_'+rcp+'.nc')['fwi']
FWI = FWI.sel(time=slice('1970-01-01','2098-12-31'))

# Linear trend function
def linear_trend(x):
    pf = np.polyfit(np.arange(1970.,2099.), x, 1)
    return xr.DataArray(pf[0])
# p-values function
def pvalue(x):
    r, p = stats.pearsonr(np.arange(1970.,2099.), x)
    return xr.DataArray(p)

# Linear trend of the number of days with FWI above the threshold
trend = xr.apply_ufunc(linear_trend,  # function to apply
                         FWI,
                         input_core_dims=(('time',),),
                         vectorize=True) 
trend.to_netcdf(path_outs+'FWI_ndays_trend_'+rcp+'.nc')

# p-values of the number of days with FWI above the threshold
pvalues = xr.apply_ufunc(pvalue,  # function to apply
                         FWI,
                         input_core_dims=(('time',),),
                         vectorize=True) 
pvalues.to_netcdf(path_outs+'FWI_ndays_pvalues_'+rcp+'.nc')